      SUBROUTINE KURV1 ( N, X, Y, SLP1, SLPN, XP, YP, TEMP, S,
     $ SIGMA )
C*
C*              *******************************************
C*              *                                         *
C*              *                   KURV1                 *
C*              *                                         *
C*              *******************************************
C*
C*      KURV1 - CLINE'S ROUTINE TO SET UP PIECEWISE CUBIC SPLINE UNDER
C*           TENSION, FOR EVALUATION BY ROUTINE KURV2.
C*
C*      AUTHOR - A. K. CLINE
C*               NATIONAL CENTER FOR ATMOSPHERIC RESEARCH
C*               PO BOX 1470
C*               BOULDER, COLORADO     80302
C*
C*      CODED BY - A. E. RAGOSTA        415-694-6235
C*               TR18
C*               AMES RSCH CTR
C*               MOFFETT FIELD, CALIF  94035
C*
C*      INPUT ARGUMENTS
C*           N     - NUMBER OF POINTS TO BE FIT
C*           X     - ARRAY OF X VALUES
C*           Y     - ARRAY OF Y VALUES
C*           SLP1  - SLOPE AT FIRST POINT (DEGREES, COUNTER-CLOCKWISE FROM
C*                       POSITIVE X AXIS)
C*           SLPN  - SLOPE AT LAST POINT (DEGREES, ...)
C*           SIGMA - TENSION FACTOR (IF THIS VALUE IS NEGATIVE, THE END
C*                     POINT SLOPES WILL BE CALCULATED, IF POSITIVE, THEY
C*                     SHOULD BE INPUT IN SLP1 AND SLP2.  A TYPICAL VALUE
C*                     IS 1. )
C*           TEMP  - SCRATCH WORK AREA
C*
C*      OUTPUT ARGUMENTS
C*           XP    - CURVATURE PARAMETERS FOR KURV2
C*           YP    - CURVATURE PARAMETERS FOR KURV2
C*           S     - ARC LENGTH OF CURVE
C*
C*      COMMON BLOCKS REFERENCED :
C*           NONE
C*
C*      FILES REFERENCED :
C*           NONE
C*
C*     EXTERNAL SUBPROGRAMS REFERENCED :
C*          SIN, SQRT, COS, ABS, FLOAT, EXP, ATAN2
C*
C*      VERSION I        29 DEC 1981
C*
C***********************************************************************
C*
      DIMENSION X(N), Y(N), YP(N), XP(N), TEMP(N)
C
      DEGRAD = 0.01745329
      NM1   = N - 1
      NP1   = N + 1
      DELX1 = X(2) - X(1)
      DELY1 = Y(2) - Y(1)
      DELS1 = SQRT ( DELX1**2 + DELY1**2 )
      DX1   = DELX1 / DELS1
      DY1   = DELY1 / DELS1
C
C --- CALCULATE SLOPES IF REQUESTED
C
      IF ( SIGMA .LT. 0. ) GO TO 70
      SLPP1 = SLP1 * DEGRAD
      SLPPN = SLPN * DEGRAD
   10 XP(1) = DX1 - COS ( SLPP1 )
C
C --- SET UP RIGHT HAND SIDE OF TRIDIAG
C
      YP(1) = DY1 - SIN ( SLPP1 )
      TEMP(1) = DELS1
      S     = DELS1
      IF ( N .EQ. 2 ) GO TO 30
      DO 20 I = 2, NM1
         DELX2 = X(I+1) - X(I)
         DELY2 = Y(I+1) - Y(I)
         DELS2 = SQRT ( DELX2**2 + DELY2**2 )
         DX2   = DELX2 / DELS2
         DY2   = DELY2 / DELS2
         XP(I) = DX2 - DX1
         YP(I) = DY2 - DY1
         TEMP(I) = DELS2
         DELX1 = DELX2
         DELY1 = DELY2
         DELS1 = DELS2
         DX1   = DX2
         DY1   = DY2
         S     = S + DELS1
   20    CONTINUE
   30 XP(N)  = COS ( SLPPN ) - DX1
      YP(N)  = SIN ( SLPPN ) - DY1
C
C --- DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS ( SIGMA ) * FLOAT ( N-1 ) / S
C
C --- FORWARD ELIMINATION ON TRIDIAGONAL
C
      DELS   = SIGMAP * TEMP(1)
      EXPS   = EXP ( DELS )
      SINHS  = .5 * ( EXPS - 1./EXPS )
      SINHIN = 1./( TEMP(1) * SINHS )
      DIAG1  = SINHIN * ( DELS * .5 * ( EXPS + 1./EXPS ) - SINHS )
      DIAGIN = 1./DIAG1
      XP(1)  = DIAGIN * XP(1)
      YP(1)  = DIAGIN * YP(1)
      SPDIAG = SINHIN * ( SINHS - DELS )
      TEMP(1) = DIAGIN * SPDIAG
      IF ( N .EQ. 2 ) GO TO 50
      DO 40 I = 2, NM1
         DELS   = SIGMAP * TEMP(I)
         EXPS   = EXP ( DELS )
         SINHS  = .5 * ( EXPS - 1./EXPS )
         SINHIN = 1./( TEMP(I) * SINHS )
         DIAG2  = SINHIN * ( DELS * ( .5 * ( EXPS + 1./EXPS )) - SINHS )
         DIAGIN = 1./( DIAG1 + DIAG2 - SPDIAG * TEMP(I-1))
         XP(I)  = DIAGIN * ( XP(I) - SPDIAG * XP(I-1))
         YP(I)  = DIAGIN * ( YP(I) - SPDIAG * YP(I-1))
         SPDIAG = SINHIN * ( SINHS - DELS )
         TEMP(I) = DIAGIN * SPDIAG
         DIAG1  = DIAG2
   40    CONTINUE
   50 DIAGIN = 1./( DIAG1 - SPDIAG * TEMP(NM1))
      XP(N)  = DIAGIN * ( XP(N) - SPDIAG * XP(NM1))
      YP(N)  = DIAGIN * ( YP(N) - SPDIAG * YP(NM1))
C
C --- PERFORM SUBSTITUTIONS FOR COEFFICIENTS
C
      DO 60 I = 2, N
         IBAK = NP1 - I
         XP(IBAK) = XP(IBAK) - TEMP(IBAK) * XP(IBAK+1)
         YP(IBAK) = YP(IBAK) - TEMP(IBAK) * YP(IBAK+1)
   60    CONTINUE
      GO TO 90
C
C --- SECOND ORDER INTERPOLATION FOR ENDPOINTS
C      (IF NO SLOPES SPECIFIED)
C
   70 IF ( N .EQ. 2 ) GO TO 80
      DELS2 = SQRT (( X(3) - X(2))**2 + ( Y(3) - Y(2))**2 )
      DELS12 = DELS1 + DELS2
      C1    = -(DELS12 + DELS1)/DELS12/DELS1
      C2    = DELS12/DELS1/DELS2
      C3    = -DELS1 / ( DELS12 * DELS2 )
      SX    = C1 * X(1) + C2 * X(2)  +  C3 * X(3)
      SY    = C1 * Y(1)  +  C2 * Y(2)  +  C3 * Y(3)
      SLPP1 = ATAN2 ( SY, SX )
      SLP1  = SLPP1 / DEGRAD +180.
      DELNM1= SQRT (( X(N-2) - X(NM1))**2 + ( Y(N-2) - Y(NM1))**2 )
      DELN  = SQRT (( X(NM1) - X(N))**2 + ( Y(NM1) - Y(N))**2 )
      DELNN = DELNM1 + DELN
      C1    = ( DELNN + DELN ) / ( DELNN * DELN )
      C2    = -DELNN / ( DELN * DELNM1 )
      C3    = DELN / ( DELNN * DELNM1 )
      SX    = C3 * X(N-2)  +  C2 * X(NM1)  +  C1 * X(N)
      SY    = C3 * Y(N-2)  +  C2 * Y(NM1)  +  C1 * Y(N)
      SLPPN = ATAN2 ( SY, SX )
      SLPN  = SLPPN / DEGRAD
      IF ( SLPN .LT. 0. )SLPN = SLPN + 360.
      GO TO 10
C
C --- TWO POINTS ONLY, RETURN A STRAIGHT LINE
C
   80 XP(1) = 0.
      XP(2) = 0.
      YP(1) = 0.
      YP(2) = 0.
      SLP1 = ATAN2 ((Y(2)-Y(1)),(X(2)-X(1))) / DEGRAD
      SLPN = SLP1
      IF ( SLPN .LT. 0. ) SLPN = SLPN + 360.
      SLP1 = SLP1 + 180.
   90 RETURN
      END
C
C---END KURV1
C
